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Abstract:  An  algorithm  is  presented  for  computing  a  column  permutation  II  and  a  QR- 
factorization  ATI  =  QR  of  an  m  by  n  (m  >  n)  matrix  A  such  that  a  possible  rank  deficiency 
of  A  will  be  revealed  in  the  triangular  factor  R  having  a  small  lower  right  block.  For  low  rank 
deficient  matrices,  the  algorithm  is  guaranteed  to  reveal  the  rank  of  A  and  the  cost  is  only  slightly 
more  than  the  cost  of  one  regular  QR-factorization.  A  posteriori  upper  and  lower  bounds  on  the 
singular  values  of  A  are  derived  and  can  be  used  to  infer  the  numerical  rank  of  A. 
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1.  Introduction 

A  very  useful  factorization  of  an  m  by  n  (m  >  n)  matrix  A  is  the  QR-factorization,  given 
by  All  =  QR,  where  n  €  f?n*n  is  a  permutation  matrix,  Q  €  Rmxn  has  orthogonal  columns  and 
satisfies  QTQ  =  /„  and  R  €  Rnxn  is  upper  triangular. 

If  A  has  full  rank,  then  R  is  nonsingular.  In  many  applications  in  which  A  is  nearly  rank 
deficient,  it  is  desirable  to  select  the  permutation  II  so  that  the  rank  deficiency  is  exhibited  in  R 
having  a  small  lower  right  block  [10].  For  if  R  is  partitioned  as 


where  R22  is  r  by  r,  then  it  is  easy  to  show  that  <7„-r+i(A)  <  ||f?22||2>  where  we  have  used  the 
notation  <r,-  to  denote  the  »-th  singular  value  of  A,  with  <ri  >  02  >  ...  >  cr„.  Therefore,  if  1 1 /?22 1 1 2 
is  small,  then  A  has  at  least  r  small  singular  values  and  thus  is  close  to  being  rank  r  deficient.  The 
converse,  unfortunately,  is  not  true.  In  other  words,  if  A  has  r  small  singular  values,  then  it  is  not 
guaranteed  that  a  given  QR-factorization  of  A  has  a  small  ||i?22||2,  as  the  following  example  shows. 

Example  1  :  Let  A„  €  Rnxn  be  the  matrix  of  order  n  illustrated  below  for  n  =  4  : 

<1  — c  -c 

,01  — c 

An  =  dia$(l,s,s2,...,srt  )  .  . 

^0  ...  0 

where  s  and  c  satisfy  s2  +  c2  =  1.  For  n  =  50, c  =  .2,  <r„(An)  »  10-4.  On  the  other  hand,  A„  is  its 
own  QR-factorization  and  obviously  has  no  small  R22  block  for  any  value  of  r. 

Besides  being  able  to  reveal  rank  deficiency  of  A,  a  QR-factorization  with  a  small  R2 2  block 
is  very  useful  in  many  applications,  such  as  in  rank  deficient  least  squares  computation  [7]  and 
in  the  subset  selection  problem  [6,  7].  Therefore  a  variety  of  techniques  have  been  proposed  to 
compute  it.  Since  the  QR  factorization  is  essentially  unique  once  the  permutation  II  is  fixed,  these 
techniques  all  amount  to  finding  an  appropriate  column  permutation  of  A.  Perhaps  the  most  well- 
known  of  these  is  the  column  pivoting  strategy  [2].  Although  this  strategy  is  usually  very  effective 
in  producing  a  triangular  factor  R  with  a  small  ||i?22||>  very  little  is  known  in  theory  about  its 
behaviour  and  it  can  fail  on  some  matrices  [5].  For  example,  column  pivoting  leaves  the  matrices 
An  in  Example  1  unchanged  and  fails  to  produce  a  small  ||i?22||-  While  such  examples  are  rare  in 
practice,  it  is  still  of  fundamental  interest  to  understand  what  distinguishes  column  permutations 
of  A  that  produce  rank-revealing  QR-factorizations  from  permutations  that  don’t,  and  to  construct 
algorithms  that  can  identify  them.  In  this  paper,  we  present  an  algorithm  that  is  guaranteed  to 
work  for  low  rank  deficient  matrices  and  costs  only  slightly  more  than  one  regular  QR-factorization 
of  A.  In  addition,  empirical  evidence  shows  that  it  works  for  most  higher  rank  deficient  matrices 
as  well.  It  is  important  to  note  that  even  if  the  singular  value  decomposition  (SVD)  of  A  is  known, 
it  is  still  not  obvious  how  to  compute  such  a  rank  revealing  QR-factorization.  In  [6],  an  algorithm 
is  presented  for  solving  the  subset  selection  problem  and  is  based  on  first  computing  the  SVD  of 
A.  Our  algorithm  does  not  require  computing  the  SVD  of  A.  In  the  rank  one  deficient  case,  the 
two  algorithms  produce  the  same  results  but  for  higher  dimensional  rank  deficient  matrices,  they 
are  different. 

In  Section  2,  we  shall  first  study  the  rank  one  deficiency  case.  In  Section  3,  we  present 
an  algorithm  for  the  general  rank  deficiency  case  and  derive  a  priori  bounds  on  the  size  of  the 


computed  ||i?22||  and  a  posteriori  upper  and  lower  bounds  on  the  singular  values  of  A.  In  Section 
4,  we  discuss  how  this  algorithm  can  be  implemented  efficiently.  Some  numerical  examples  will  be 
given  in  Section  5. 

We  shall  use  the  notation  x,-  or  (a?);  to  denote  the  i-th  component  of  a  vector  x,  and  a, 7  to 
denote  the  (*,  j)— th  element  of  matrix  A. 

2.  Revealing  Rank  One  Deficiency 

Assume  that  A  is  nearly  rank  one  deficient.  We  would  like  to  find  a  column  permutation 
of  A  such  that  the  resulting  QR- factorization  has  a  small  pivotal  element  rnn.  It  turns  out  that 
this  permutation  can  be  found  by  inspecting  the  size  of  the  elements  of  the  singular  vector  of  A 
corresponding  to  the  smallest  singular  value  crn.  This  procedure  was  first  pointed  out  in  [6]. 

Theorem  2.1. 

Suppose  that  we  have  a  vector  x  €  Rn  with  ||x||2  =  1  such  that  ||Ax||2  =  e  and  let  II  be  a 
permutation  such  that  ifUTx  =  y,  then  |y„|  =  j  |y 1 1^ .  Then  if  All  =  QR  is  the  QR- factorization  of 
All,  then 

knn|  <  \fnt. 


Proof. 

The  proof  of  this  theorem  was  stated  as  Exercise  P6.4-4  in  [7].  We  shall  prove  it  here.  First 
we  note  that  since  |y„|  =  Hj/H^  and  ||y||2  =  ||x||2  =  1,  we  have  |y„|  >  Next  we  have 

QtAx  =  QTAfmTx  =Ry=( 

\rnnyn 

Therefore, 

€  =  ||Ax||2  =  ||<?rAx||2  =  ||i?y||2  >  |r„„y„|, 
from  which  the  result  follows. 


Let  v  €  Rn  with  ||v„||  =  1  be  the  right  singular  vector  of  A  corresponding  to  the  smallest 
singular  value  <r„.  Then  we  have 

II  Av||2  =  crn. 

Therefore,  by  Theorem  2.1,  if  we  define  the  permutation  II  by 

|(nrv)L  =  IML 

then  All  has  a  QR-factorization  with  a  pivot  rnn  at  least  as  small  as  \fnon  in  absolute  value. 

Since  only  the  permutation  II  is  needed,  it  is  not  necessary  to  compute  the  SVD  of  A  in  order 
to  find  v  exactly.  In  practice,  one  can  use  a  few  steps  of  inverse  iteration  [3,  9]  to  compute  an 
approximation  to  v  from  which  the  permutation  IT  can  be  determined.  In  the  more  interesting  case 
where  <r„  <C  the  inverse  iteration  should  converge  rapidly.  This  suggests  a  two-pass  algorithm 
in  which  one  first  computes  any  QR-factorization  of  A,  then  performs  inverse  iteration  with  R  to 
find  an  approximate  v,  determines  II  and  then  computes  the  QR  factorization  of  All.  The  above 
procedure  is  similar  to  a  two-pass  algorithm  derived  in  [4]  for  computing  a  LU  factorization  of  a 
square  matrix  B  with  a  small  n-th  pivot.  In  fact,  the  two  algorithms  are  very  much  related  but 
we  will  not  pursue  that  here. 


3.  Revealing  Higher  Dimensional  Rank  Deficiency 

In  this  section,  we  consider  the  case  where  A  is  nearly  rank  r  deficient,  with  r  >  1.  Our  goal 
is  to  find  a  permutation  II  such  that  if 

AH  =  QK  j£) 

is  the  QR-factorization  of  >111,  with  P22  €  Rrxr  then  ||J?22 1|  is  small  in  some  norm. 

A  natural  way  to  extend  the  one  dimensional  result  of  Section  2  is  to  repeatedly  apply  the  one 
dimensional  algorithm,  for  r  =  1, 2, . . . ,  to  i?n,  the  leading  principal  triangular  part  of  R.  Suppose 
that  we  have  already  isolated  a  small  r  by  r  block  i?22-  To  isolate  a  small  (r  +  1)  by  (r  +  1)  block, 
we  compute^  using  the  one  dimensional  algorithm  given  in  Section  2,  a  permutation  P  such  that 
f?nP  =  Qii?n  is  the  QR-factorization  of  RnP  and  where  the  (n  —  r, n  —  r)-th  element  of  Rn  is 
small.  Then  with 

fisn(o  l) 

5) 

it  can  be  easily  verified  that 

'iff=o(so1 

is  the  QR-factorization  of  All. 

The  above  procedure  can  be  summarized  in  the  following  algorithm  : 

Algorithm  RRQR(r) 

Purpose:  Computes  a  permutation  II  and  a  QR-factorization 

An  =  <3*  =  <?(*'  *;) 

with  a  small  f?22  €  Rrxr  to  reveal  a  possible  rank  r  deficiency  in  A. 


Compute  any  QR-factorization  of  A  :  All  =  QR. 
Initialize  W  €  Rnxr  to  zero. 

For  i  =  n,  n  -  I, . . . ,  n  —  r  -f  1,  do 


1.  f?n  *—  leading  t  x  i  block  of  R. 

2.  Compute  the  singular  vector  v  €  R'  of  R\\  corresponding  to  <rmin(Pii)  with  1 1 v| I2  =  1  and  set 
bi  =  ^min(-^ll)* 

3.  Compute  a  permutation  P  €  R,x>  such  that  |(PTv),|  =  | IP*7, vj  1^ - 


4.  Assign  v  =  G  Rn  to  the  i-th  column  of  IV. 

5.  Compute  W  <—  PTW,  where  P=  ^  . 


6.  Compute  the  QR-factorization  :  J?nP  =  Qif?u. 
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*■<>-<>(%' 

End  For 

The  matrix  W  is  used  to  store  the  singular  vectors  corresponding  to  the  smallest  singular  values 
of  the  successive  leading  triangular  blocks  i?n-  It  is  updated  in  the  above  algorithm  for  theoretical 
reasons  only  (see  Lemma  3.2)  and  need  not  be  explicitly  incorporated  in  an  actual  implementation. 

For  the  above  algorithm  to  produce  the  desired  QR- factorization,  we  must  prove  the  following 
two  assertions  : 

1.  At  step  2,  J?n  has  a  small  singular  value  so  that  the  (i,  t)-th  element  of  /?n  is  guaranteed  to 
be  small. 

2.  At  step  9,  the  last  (i.e.,  the  (n  -  i)-th)  row  of  Q\ Rn  is  small. 

If  these  two  assertions  are  true,  then  the  lower  (n  —  i  +  1)  by  (n  —  i  +  1)  lower  block  of  R  in  step  9 
is  small  and  we  have  the  desired  QR-factorization. 

We  shall  next  prove  more  precise  versions  of  these  two  assertions.  With  regard  to  the  first 
assertion,  we  have  the  following  lemma. 

Lemma  3.1.  Let  B  €  Rnxk  be  a  matrix  containing  any  subset  of  k  columns  of  A,  then 


^min(fl)  —  &k(B)  ^  ^k(A). 


Proof.  Follows  from  the  variational  characterization  of  singular  values.  See  [7], 


Since  Ru  is  the  first  i  columns  of  QTAU,  by  the  Lemma  3.1,  we  have 

^minCRll)  ^  &i(Q  All)  =  (Tj (A). 

This  guarantees  that  f?n  has  a  small  singular  value  if  ct,(A)  is  small,  i.e.,  if  the  algorithm  has  not 
yet  “captured”  the  complete  approximate  null  space  of  A. 

Proving  the  second  assertion  is  more  subtle.  We  shall  first  need  to  establish  some  properties 
of  the  matrix  W  employed  in  the  algorithm. 

Lemma  3.2.  The  matrix  W  =  [w„_r+i,  . . . ,  w„]  €  Rnxr  computed  by  Algorithm  RRQR(r)  satisfies 
the  following  properties  : 

For  i  =  n,  n  -  1, . . . ,  n  —  r  +  1 

1.  |H|,  =  i. 

2.  (wi)j  =  0  for  j  >  i. 

3.  IK), |  =  IKIL  >  ^7- 

4.  || AITw,||2  =  h  <  <MA). 

Proof.  We  shall  prove  the  lemma  by  induction  on  r.  The  lemma  is  obviously  true  for  the  case  r  =  1 
because  of  the  one  dimensional  result  in  Theorem  2.1.  Now  assume  that  the  lemma  is  true  for  r  =  k. 


That  is,  after  k  steps  of  Algorithm  RRQR(fc),  the  vectors  wj,  for  j  =  n,  n  —  1, . . . ,  n  —  k+ 1,  satisfies 
the  four  properties  stated  in  the  lemma.  We  shall  show  that  the  vectors  Wj,  for  j  =  n, . . .  ,n  -  k, 
computed  by  Algorithm  RRQR(fc  +  1)  also  satisfy  these  properties.  First  observe  that  Algorithm 
RRQR(fc+ 1)  is  simply  Algorithm  RRQR(lc)  followed  by  one  more  pass  through  the  loop  with  index 
i  =  n  —  k.  At  step  5  of  this  last  pass,  the  matrix  W  is  updated  by 

W  <-PT{v,wn-k+l,...,wn) 

={PTv,  PTwn.k+i PTwn), 

where  by  the  induction  hypothesis,  the  vectors  Wj,  j  =  n,. . .  ,n  -  k  +  1  satisfies  the  properties 
stated  in  the  lemma.  Note  that  since  PT  is  a  permutation  affecting  only  the  first  n  —  k  components 
of  each  of  the  wj  vectors,  the  updated  vectors  wj,  for  j  =  n, ...,»  —  ifc  +  1,  automatically  satisfy  the 
first  three  properties.  Moreover,  v  is  chosen  so  that  wn-k  =  PTv  also  satisfy  these  three  properties. 
To  verify  the  fourth  property  for  j  =  n, . . . ,  n  -  k  +  1,  note  that  the  permutation  for  Algorithm 
RRQR(&  +  1)  is  IIP  and 

llAnP^Ha  =  ||AIIPPr%||2  =  ||AIIuv||2  <  <r> 
by  the  induction  hypothesis.  Finally,  for  j  =  n  -  k,  we  have 
||AIIP«/„-*||2  =  ||AnPPrv||2  =  ||AIIu||2 

<Ro  SOGOI, 

=  II-Riiv||2 

==  ^minC^ll) 

=  f>n-k 

<  <rn-k(A) 

where  the  last  bound  is  arrived  at  by  using  Lemma  3.1,  and  the  fact  that  R  has  the  same  singular 
values  as  A.  This  completes  the  inductive  proof. 

I 

We  are  now  ready  to  state  our  main  results. 

Theorem  3.1.  Let  the  matrix  W  €  Rnxr  computed  by  Algorithm  RRQR(r)  be  partitioned  as 
WT  =  (W1r,  Wj),  where  W2  €  Rrxr  is  upper  triangular  and  nonsingular.  Then  the  QR-factorization 
All  =  QR  as  computed  by  Algorithm  RRQR(r)  satisfies 

||P22||2<<Tn_r+l||^2-l||2. 

Proof.  Denote  the  columns  of  W  by  [w„_r+i,  - Define  y  =  <5T.4IIu;„_r+i  e  Rnxr.  From 
property  4  of  Lemma  3.2,  we  have 

||y||2  =  S„-r+ 1  <  (Tn-r+i  ■ 

On  the  other  hand,  with  e\  =  (1,0,  ■  •  •  ,0)r  e  Rr,  we  have 

»  =  QTAmVe,  =  BH'„  =  (  R"W£J"W')  eu 
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from  which  follows  that 


Combining  the  above,  we  get 


ll!/||2  >  11-^22^2^1 1|2  > 


1 1  -^22 1 1 2 


&n—r+ 1  —  ^n-r+1  ^ 


II-R22II2 


from  which  the  desired  result  follows. 


This  theorem  states  that  HR22II2  is  bounded  by  the  r-th  smallest  singular  values  of  A,  provided 
that  not  to°  iar8e-  From  Lemma  3.2,  we  have  that  W2  is  upper  triangular,  each  one  of  its 

columns  has  unit  Euclidean  length  and  is  obviously  nonsingular.  Moreover,  Algorithm  RRQR  can 
be  interpreted  as  trying  to  produce  a  well-conditioned  Wi  by  putting  the  largest  element  in  absolute 
value  of  each  column  on  the  diagonal.  In  this  sense,  it  is  analogous  to  the  strategy  proposed  in 
[6],  in  which  the  permutation  of  II  is  obtained  via  a  QR-factorization  with  column  pivoting  applied 
to  isolate  a  well-conditioned  subset  of  the  rows  of  the  matrix  formed  by  the  r  singular  vectors 
corresponding  to  the  r  smallest  singular  values  of  A.  The  present  algorithm  has  the  advantage  that 
the  SVD  of  A  need  not  be  computed  and  that  an  a  priori  bound  for  ||i?22||  can  be  derived  from 
an  a  priori  upper  bound  for  ||W'£"1||2.  In  fact,  an  a  priori  bound  for  individual  elements  of  R22  is 
possible. 

Theorem  3.2.  Algorithm  RRQR(r)  computes  a  permutation  II  and  a  QR-factorization  of  A  given 
by  All  =  QR  where  the  elements  of  the  lower  r  by  r  upper  triangular  block  of  R  satisfy  : 

i- 1 

k—i 

<  2 J~'<7jy/n  for  n  —  r  <  i  <  j  <  n.  (3.2) 

Proof.  By  employing  Lemma  3.2,  we  have,  for  n  —  r  <  i  <  j  <  n, 

j 

'52rik{'»j)k  • 

k=i 

Isolating  the  k  =  j  term  in  the  sum  and  rearranging  terms,  we  get 


aj  >  ||AIItnj2  =  WRuJjh  >  |(i?«v).|  = 


(3.1) 


>-i 

M«V)y|  <  Oj  +  ^2 \rik{wj)k\  ■ 
k=i 

Since  by  Lemma  3.2,  |(w,-),-|  =  ||(w.)||00  <  -7;,  we  have 

V> 

1-1 

ik=i 

Solving  this  recurrence  in  the  index  j,  we  get  the  bound  given  in  (3.2).  Using  the  bound  Ok\/k  < 
<Tl  V/n  in  each  term  in  the  sum  in  (3.2),  we  get  the  bound  in  (3.2) 


The  term  2J~'  in  the  bound  (3.2)  in  Theorem  3.2  indicates  that  the  bound  for  an  element  of 
f?22  is  larger  the  further  the  element  is  from  the  main  diagonal.  In  fact,  the  bound  for  the  element 
|rn_r+i,n|  grows  like  2r.  For  large  values  of  r,  this  bound  can  be  quite  large  and  overly  conservative 
for  most  problems  encountered  in  practice.  On  the  other  hand,  many  problems  in  application  have 
small  values  of  r,  and  for  these  problems,  the  bounds  in  Theorem  3.2  should  be  quite  reasonable. 
For  example,  if  r  =  4,  then  we  have,  for  n  -  4  <  t  <  j  <  n, 

max  |r,j  |  <  Wy/ncr,-. 

0 

Therefore,  any  small  singular  values  of  A  will  be  revealed  in  the  triangular  factor  R.  In  other 
words,  for  low  rank  deficient  matrices,  Algorithm  RRQR(r)  is  guaranteed  to  produce  rank  revealing 
QR-factorizations. 

4.  Bounding  the  Singular  Values 

Very  often,  it  is  desirable  to  be  able  to  infer  the  rank  of  A  from  its  QR-factorization  by 
estimating  the  small  singular  values  of  A  from  the  triangular  factor  R  [l,  8].  From  Theorem  3.1, 
we  can  easily  derive  the  following  a  posteriori  bounds  on  the  singular  values  of  A  in  terms  of  the 
matrices  R  and  W  computed  by  Algorithm  RRQR(r). 

Corollary  4.1.  Let  W  €  Rnxr  and  R  €  Rnxn  be  as  computed  by  Algorithm  RRQR(r).  Let  R}22 
and  denote  the  lower  right  j  by  j  upper  triangular  blocks  of  R  and  W  respectively.  Then,  for 
1  <  j  <  r, 

||(VF',)~1||  ~  ^ n-/+ 1  -  en-j+l  <  Il-^22ll2  —  ^fW+lIK^)  1  lb- 

The  quantities  and  H-R^lh  are  computable  a  posteriori  lower  and  upper  bounds  respec¬ 

tively  for  the  singular  value  <rn_y+i.  If  ||(IV^)_1||2  is  not  much  larger  than  the  optimal  value  of 
one,  then  the  above  Corollary  states  that  these  bounds  should  be  very  tight  and  one  can  therefore 
infer  the  numerical  rank  of  A.  In  practice,  the  value  of  HR22II2  can  be  estimated  by  more  easily 
computable  norms  of  R22,  such  as  the  1-norm,  the  00— norm  and  the  Frobenius  norm.  These  a 
posteriori  bounds  also  enable  Algorithm  RRQR  to  be  implemented  in  an  adaptive  manner:  the 
algorithm  can  be  terminated  if  the  lower  bounds  indicate  that  all  the  small  singular  values  of  A 
have  been  revealed. 

5.  Implementation  Note 

The  main  work  of  Algorithm  RRQR  consists  of  three  parts  :  the  computation  of  the  initial  QR- 
factorization,  the  computation  of  the  singular  vector  v  of  Ru  by  inverse  iteration  at  each  iteration 
and  the  QR-factoring  of  RnP  at  each  iteration.  Ignoring  lower  order  terms,  the  first  two  takes 
n2(m  -  3)  (assuming  that  Q  needs  not  be  accumulated)  and  In2r  flops  respectively,  where  I  is  the 
number  of  inverse  iterations  used  at  each  each  iteration.  Usually,  I  =  2  is  sufficient  in  practice, 
especially  if  an  efficient  condition  number  estimator  is  used  to  select  the  initial  vector. 

The  re-factoring  of  RnP ,  if  done  naively,  say  by  using  Householder  transformations,  would  be 
extremely  inefficient.  Due  to  the  special  structure  of  the  problem,  a  much  more  efficient  procedure 
can  be  derived  by  using  Given’s  rotation  instead.  The  column  of  Rn  to  be  permuted  to  the  last 
column  can  be  swapped  with  adjacent  columns  one  at  a  time  until  the  last  column  is  reached.  For 
example,  if  column  1  is  to  be  permuted  to  column  k,  we  can  swap  the  pairs:  (1,2),  (2,3), . . . ,  (fc— I,  it) 
successively.  With  each  swap,  exactly  one  nonzero  is  introduced  in  the  lower  triangular  part  of  f?u, 


Page  8 


which  can  then  be  annihilated  by  applying  an  appropriate  Given’s  rotation  from  the  left.  In  most 
applications,  these  extra  Given’s  rotation  can  be  stored  in  factored  form  and  do  not  have  to  be 
applied  to  the  orthogonal  factor  Q  of  the  initial  QR-factorization.  Assuming  the  worst  case  in  each 
iteration,  i.e.,  that  the  first  column  of  Ru  is  to  be  permuted  to  the  last  every  time,  this  takes  2 n2r 
flops.  Therefore,  the  total  work  for  Algorithm  RRQR(r)  is  given  by  : 

W  (r)  =  n2(m  —  x)  +  +  2n2r 

=  n2(m  —  t)  +  4n2r,  assuming  that  1  =  2. 

O 

Therefore,  if  r  4C  n,  the  extra  work  performed  after  the  initial  QR  factorization  is  of  a  lower  order 
and  negligible,  especially  if  m  »  n.  Even  in  the  extreme  case  of  r  =  n,  W{n)  =  mn 2  +  ^-n3,  which 
is  still  smaller  than  the  cost  of  2mn2  +  4n3  for  computing  the  SVD. 

It  is  interesting  to  note  that  the  overhead  for  the  column  pivoting  strategy  is  O(mn)  while  the 
overhead  for  RRQR  is  0(n2r).  Therefore,  if  r  is  small  and  m  >>  n,  the  overhead  of  RRQR  is  less 
than  that  for  column  pivoting. 


6.  Numerical  Examples 

In  this  section,  we  present  a  few  numerical  examples  to  illustrate  the  theory  developed  in 
the  earlier  sections.  All  computations  are  done  on  a  VAX/780  in  single  precision,  with  a  relative 
machine  precision  of  about  10-7.  The  initial  QR-factorization  is  performed  with  no  column  per¬ 
mutation  and  the  inverse  iterations  are  performed  with  1=5  with  the  initial  vector  chosen  to  be 
(1,0, 1,0,  •  •  • ,  1,0)T.  For  each  of  the  examples,  we  compute  the  QR-factorization  with  Algorithm 
RRQR,  compute  the  lower  and  upper  a  posteriori  bounds  and  H.R22II2  given  in  Corollary 

4.1,  and  compare  these  to  the  true  singular  values  <rn~j+ 1  a-nd  to  the  upper  bounds  obtained  by 
the  regular  QR-factorization  with  and  without  column  pivoting. 

The  first  example  is  the  matrix  A50  given  earlier  in  the  introduction.  It  turns  out  that,  due  to 
round-off  errors,  QR  with  column  pivoting  actually  pivots  with  this  matrix.  We  therefore  use  the 
following  modified  matrix  instead 


An  =  A„  +  diag(ne,  ( n  -  l)e,-  •  ,e), 


where  e  =  10-6.  This  small  perturbation  only  change  the  singular  values  of  the  matrix  by  a  small 
relative  amount  and  forces  no  column  pivoting.  The  results  are  shown  in  Table  1,  in  which  the 
column  permutations  produced  by  RRQR  and  QR  with  column  pivoting  are  also  presented.  We 
see  that  in  this  case,  Algorithm  RRQR  succeeds  in  revealing  the  one  dimensional  null  space  of 
the  matrix  whereas  QR  with  column  pivoting  fails  to  do  so.  Since  the  values  of  ||(W^)-1||2  are 
very  close  to  one,  the  a  posteriori  bounds  are  very  tight  and  there  is  no  trouble  in  identifying 
the  numerical  rank.  In  fact,  since  the  lower  bound  for  U49  is  clearly  not  small,  one  can  confidently 
truncate  the  algorithm  after  only  two  iterations.  Note  that  since  the  lower  bounds  <5,  ’s  are  computed 
only  approximately  by  inverse  iteration,  they  can  actually  be  slightly  larger  than  the  true  singular 
values.  But  this  slight  inacurracy  does  not  affect  the  ability  to  infer  the  rank. 

The  second  example  is  the  class  of  matrices  defined  as  follows.  Let  H„  =  I  —  j^eer,  where 
e  =  (1,1,...,  l)r.  Define 

C  =  H™(o)H'0' 

where  D  is  a  diagonal  matrix.  Since  H„  is  orthogonal,  the  singular  values  of  C  are  the  diagonal 
entries  of  D.  The  results  for  two  different  D' s  with  well-defined  numerical  rank  are  presented  in 
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i 

RRQR 

lower 

true 

<7,' 

RRQR 

upper 

K 

pivoted  QR 
upper 

K 

regular  QR 
upper 

50 

9.29E-05 

0.0001 

0.0002 

1 

0.3678 

50 

0.3678 

49 

0.4349 

0.4112 

0.4538 

48 

0.4112 

49 

0.4112 

48 

0.4262 

0.4222 

0.4937 

49 

0.4227 

48 

0.4227 

47 

0.4531 

0.4327 

0.4967 

46 

0.4351 

47 

0.4351 

46 

0.4440 

0.4431 

0.5143 

47 

0.4495 

46 

0.4495 

45 

0.4722 

0.4536 

0.5179 

44 

0.4665 

45 

0.4665 

44 

0.4627 

0.4642 

0.5357 

45 

0.4868 

44 

0.4868 

43 

0.4924 

0.4749 

0.5399 

42 

0.5106 

43 

0.5106 

42 

0.4824 

0.4858 

0.5581 

43 

0.5382 

42 

0.5382 

41 

0.5137 

0.4968 

0.5627 

40 

0.5698 

41 

0.5698 

Table  1:  Results  for  A50  ( K  =  corresponding  column  of  A 
permuted  to  i-th  column  of  All) 


Tables  2-3.  In  Table  4,  we  present  an  example  where  the  numerical  rank  is  not  as  well-defined.  All 
the  results  show  that  the  a  posteriori  bounds  producd  by  RRQR  are  tight  and  that  they  reveal  the 
numerical  rank  as  predicted.  QR  with  column  pivoting  also  succeeds  on  these  examples  although 
it  is  not  possible  to  infer  the  rank  since  a  lower  bound  on  the  singular  values  are  not  available.  The 
column  permutations  produced  are  also  different  than  those  of  RRQR.  Regular  QR-factorization 
with  no  pivoting  fails  in  all  cases. 

In  all  the  above  examples,  the  actual  values  of  11(^2  )_1||2  are  not  much  larger  than  the  optimal 
value  of  one,  and  are  much  smaller  than  that  suggested  by  the  a  priori  bounds  given  in  Theorem 
3.2.  This  is  probably  typical  for  most  problems  encountered  in  practice.  In  fact,  we  have  not 
succeeded  in  finding  an  example  for  which  the  a  priori  bounds  are  achieved.  Thus,  not  only  is 
Algorithm  RRQR  guaranteed  to  work  for  low  rank  deficient  matrices,  it  will  almost  always  work 
for  high  rank  deficient  ones  as  well. 

7.  Conclusion 

In  this  paper,  we  have  presented  a  simple  and  efficient  algorithm  for  computing  a  QR- 
factorization  that  is  designed  to  reveal  the  numerical  rank  of  a  given  matrix.  Our  theory  shows 
that  the  algorithm  is  guaranteed  to  work  for  low  rank  deficient  matrices  and  numerical  experiments 
indicate  that  it  is  also  very  likely  to  work  even  in  the  high  rank  deficient  cases.  The  a  posteriori 
bounds  for  the  singular  values,  which  are  by-products  of  the  algorithm,  can  be  used  to  infer  the 
numerical  rank  of  the  matrix,  without  explicitly  computing  the  SVD  of  the  matrix.  Thus,  this  spe¬ 
cial  QR-factorization  could  be  used  reliably  and  efficiently  in  many  matrix  computations  for  which 
a  rank-revealing  QR-factorization  is  needed,  such  as  in  rank  deficient  least  squares  computations 
and  the  subset  selection  problem. 

Acknowledgement  :  The  author  would  like  to  thank  Dr.  Robert  Schreiber  for  helpful  discussions 
on  the  development  of  Algorithm  RRQR  and  Professor  Charles  Van  Loan  for  pointing  out  Theorem 
2.1  in  the  book  [7]. 
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■ 

RRQR 

lower 

true 

RRQR 

upper 

K 

pivoted  QR 
upper 

K 

regular  QR 
upper 

10 

0.0001 

0.0001 

0.0002 

3 

0.0001 

10 

0.0001 

9 

0.0001 

0.0001 

0.0002 

6 

0.0001 

9 

0.0001 

8 

0.0001 

0.0001 

0.0002 

7 

0.0002 

8 

0.0002 

7 

0.0001 

0.0001 

0.0002 

10 

0.0002 

7 

0.0002 

6 

0.0001 

0.0001 

0.0002 

9 

0.0002 

5 

1.0000 

5 

0.4472 

1.0000 

1.0000 

8 

1.0000 

6 

1.0000 

■lECTEM 

1.0000 

1.0000 

4 

1.0000 

4 

1.0000 

Table  2:  Results  for  C  with  D  =  diag(l,  1,1, 1,1, 10  4,10  4,10  4,10  4,10  4) 

■ 

RRQR 

lower 

true 

<r, i 

RRQR 

upper 

K 

pivoted  QR 
upper 

H 

regular  QR 
upper 

10 

0.0001 

0.0001 

0.0002 

5 

0.0001 

0.0001 

9 

0.0001 

0.0001 

0.0002 

2 

0.0001 

8 

0.0002 

8 

0.0001 

0.0001 

0.0002 

8 

0.0002 

6 

0.0002 

7 

0.0001 

0.0001 

0.0002 

4 

0.0002 

4 

1.0000 

6 

0.0001 

0.0001 

0.0002 

6 

0.0002 

9 

1.0000 

5 

0.4472 

1.0000 

1.0000 

10 

1.0000 

2 

1.0000 

4 

0.6325 

1.0000 

1.0000 

7 

1.0000 

7 

1.0000 

Table  3:  Results  for  C  with  D  =  diag(l,10  4 ,1,10  4,1,10  4 , 1,10  4,1,10  4) 


RRQR 


true 


RRQR 


pivoted  QR 


regular  QR 


10 

9 

8 

7 


1 


lower 

9.8E-06 

0.0001 

0.0009 

0.0076 

0.0706 

0.4479 

0.4490 

0.6334 

0.7767 

0.8947 


upper 


upper 


9.8E-06 


1.2E-05 


1.3E-05 


0.0001 


0.0001 


0.0001 


0.0009 


0.0017 


0.0017 


0.0100 


0.0263 


0.0263 


0.1000 


0.2193 


0.2193 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


1.0000 


10 


1.0000 


K 

1 

2 

3 
10 

4 

5 
9 
8 
7 

6 


upper 

4.9E-05 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


Table  4:  Results  for  C  with  D  —  diag(10  5,10  4,10  3,10  2,10  1,1,1,1,1,1) 
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